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ABSTRACT 

One dimensional versions of dissipationless cosmological N-body simulations have been 
shown to share many qualitative behaviours of the three dimensional problem. Their 
interest lies in the fact that they can resolve a much greater range of time and length 
. scales, and admit exact numerical integration. We use such models here to study how 

O f non-linear clustering depends on initial conditions and cosmology. More specifically, 

we consider a family of models which, like the three dimensional EdS model, lead for 
power-law initial conditions to self-similar clustering characterized in the strongly non- 
linear regime by power-law behaviour of the two point correlation function. We study 
how the corresponding exponent 7 depends on the initial conditions, characterized by 
the exponent n of the power spectrum of initial fluctuations, and on a single parameter 
K controlling the rate of expansion. The space of initial conditions/cosmology divides 
very clearly into two parts: (1) a region in which 7 depends strongly on both n and 
k and where it agrees very well with a simple generalisation of the so-called stable 
clustering hypothesis in three dimensions, and (2) a region in which 7 is more or less 
independent of both the spectrum and the expansion of the universe. The boundary 
in (n, k) space dividing the "stable clustering" region from the "universal" region is 
very well approximated by a "critical" value of the predicted stable clustering expo- 
nent itself. We explain how this division of the (n, k) space can be understood as a 
simple physical criterion which might indeed be expected to control the validity of 
the stable clustering hypothesis. We compare and contrast our findings to results in 
three dimensions, and discuss in particular the light they may throw on the question 
of "universality" of non- linear clustering in this context. 

Key words: Cosmological structure formation, gravitational clustering, TV-body sim- 
ulation 



1 INTRODUCTION 

Cosmological TV body simulations are the primary instru- 
ment used to make theoretical predictions for structure for- 
mation in current models of the universe. Analytical un- 
derstanding of the non-linear regime, crucial for many non- 
trivial tests of these models, remains very poor despite the 
many phenomenological results derived in simulations. Bet- 
ter understanding of the physics of non-linear clustering 
could help to better constrain and control cosmological sim- 
ulations, which despite their ever increasing power and so- 
phistication are still subject to questions about their ac- 
curacy and resolution. In this work we use simplified one 
dimensional (ID) "toy models" of the full three dimensional 
(3D) problem to study non-linear structure formation. More 
specifically we study, using a family of such models, how the 



properties of non-linear clustering depend on both initial 
conditions and cosmology. Our main result is that, for cold 
initial conditions, we observe that the large phase space of 
initial conditions and cosmology we explore breaks up into 
two parts: on the one hand, a region in which there is an 
apparent "universality" of the non- linear clustering, i.e., in 
which it depends very weakly both on initial conditions and 
on the cosmological expansion; and on the other hand, a 
region where the non-linear clustering depends strongly on 
both the initial conditions and expansion history. In the lat- 
ter region, the exponent characterizing the two point cor- 
relations is in excellent agreement with a simple analytical 
prediction which follows from the so-called stable clustering 
hypothesis, appropriately generalized to this class of mod- 
els. Further the apparently well localized boundary dividing 
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the phase space into these two disjoint regions can in fact be 
characterized using the stable clustering prediction, in terms 
of a critical value of the corresponding exponent. The non- 
linear clustering appears to be strictly universal in the limit 
in which there is no expansion, but also very well approx- 
imated for steeper initial spectra when there is expansion. 
We compare our results with those in three dimensions, and 
discuss how our results may provide a framework for under- 
standing the much discussed apparent "universality" of halo 
profiles in this context. 

The exploration of ID models as a tool for under- 
standing the physics of non-linear structure formation 
cosmological context g oes back at least as far as 



the work of iMelottl (|l982l . 119831 ) for the case of hot 
dark matter models. Different initial conditions and vari- 
ants of the m odel have been discussed by a number of 
other authors ilYano fc Goudalll998l; iMiller fc Rouetll20ol 



Aurcll & Fanclli 2002 

' 120071 ; iGabrielli et all 120091 : 



Miller & Rouct 2006; Valageas 2006; 

Miller fc Rouetl 



Miller et al, 

2010al ). There is also an extensive literature on the statisti- 



cal mechanics of finite one dimensional se l f-grav itating sys- 
tems (see e.g. IJovce fc Worrakitpoonponl l|201ll) and refer- 
ences therein), as well as some studies which use this case to 



explor e issues in cosmological N body simulations: Binncv 



2004) uses it to probe discreteness effects, and lSchulz et al 



M) 
12) 



2012J) the issue of universality in cold collapse. 



The study reported here is a direct continuation of that 
in I Joyce fc Sicardl l|201ll ). This latter paper showed, for the 
three ID models studied in the previous literature, the very 
strong qualitative similarities in the temporal development 
of clustering to that in the 3D case. For cold cosmological 
like initial conditions clustering is hierachical and driven by 
the linear amplification of fluctuations, leading, for the case 
of power law initial conditions, to self similar evolution of 
the non-linear correlations. The exponent characterising the 
non-linear clustering was, further, found to be in very good 
agreement with that predicted by the so-called "stable clus- 
tering" hypothesis. In this paper we extend this study to a 
one parameter family of models which allows continuous in- 
terpolation between the static model and the two expanding 
models previously considered. This allows us to probe fully 
the extent of validity of stable clustering, where it breaks 
down and what happens in this case. The study reveals that 
there is a fairly abrupt switchover from the validity of stable 
clustering to an apparent universality in the clustering. 

In the next section we define the class of ID models we 
study, discussing also their relation to 3D cosmological mod- 
els. This presentation improves and simplifies that given in 
IJovce fc Sicardl l|20 11]) . For the family of scal e- free models we 
focus on, we derive, more rigorously than in I Joyce fc Sicardl 
(2011), the prediction for the behaviour of an isolated struc- 
ture and use it to derive the exponents predicted for the 
non-linear two point correlation when the stable clustering 
hypothesis is made. In the following section we discuss our 
numerical simulations, and in particular how we calibrate 
them using "exact" simulations. The next section gives our 
results, and in the last section we discuss their interpretation 
and possible relation to the three dimensional problem. 



2 A FAMILY OF ID SCALE-FREE MODELS 

The class of ID models we study are defined as follows. We 
take the equations of motion in comoving coordinates of 3D 
cosmological JV body simulations, and simply replace the 
3D Newtonian gravitational interaction with the ID gravi- 
tational interaction. This gives ID equations of motion 



HF +2H Ht 



(1) 



where Xi are the particle positions on the line. The super- 
script 'J' in the force term indicates that the sum, which 
extends over an infinite distribution of masses with non-zero 
mean mass density, is regularized, just as in 3D, by subtrac- 
tion of the contribution of the mean mass density, i.e., the 
force term is sourced only by fluctuations about the mean 
density. As we will discuss further below this sum can be 
written explicitly in a simple manner. 

The coupling g is the ID analogy of Newton's constant 
(multiplied by the particle mass). To establish a more direct 
relation between the ID and 3D models, one can consider 
the "particles" in the ID system to represent infinite parallel 
sheets of zero thickness embedded in 3D. In this case one has 
the identification g — 2irGY, where £ is the surface mass 
density of the sheets. The mean (comoving) mass density po 
of the corresponding 3D universe and the number density no 
of the ID system are then related as po = Eno, and therefore 



gn = 2nGp ■ 



(2) 



For any given 3D cosmological model which provides a scale 
factor a(t) and associated Hubble constant H(t) we then 
have a unique corresponding ID model. The case a(t) = 1 
(and H (t) — 0) on the other hand, defines a model analogous 
to a 3D universe without expansion (which is well defined 
but does not correspond to a cosmological model derived 
from general relativity). 

The one parameter family of models we will study is 
a continuous interpolation between the ID model obtained 
using the 3D EdS cosmology, and the static model. To define 
them precisely, and motivate their choice, it is convenient 
to change time variable in Eq. {TJ defining r = J -0^ . This 
gives the equations of motion in the form 



d 2 Xi . .dxi \ ^ , 

-J^+nr)— = -.9^sgn( 



(3) 



i.e., in which all effects of the cosmology appear only as a 
fluid type damping term, with 

1 3/9. „ 1 da 



T = -a 3/2 H - 

2 2a dr 



(4) 



For a 3D EdS universe with comoving mass density po, i.e., 
a = (i/io) 2 ^ 3 where to = 1/ y/QnGpo, we obtain a Y which is 
independent of time and given by 

1 



3to 



x/27rG/9o/3 = 



(•>) 



where we have used Eq. @ to derive the last equality. In 
a ACDM cosmology, on the other hand, F is approximately 
constant through the matter dominated era and increases 
slightly at late times. 

We consider here simply the set of ID models in which 
r takes any constant value, and the focus of our study is 
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the dependence of non-linear clustering on this parameter, 
as well as on the initial conditions. T can simply be con- 
sidered as a control parameter for understanding the role 
of expansion in the determination of the properties of the 
non-linear clustering. While most previous studies have con- 
sidered either the EdS model (also known as the "quintic" 
model for reasons which will be recalled below), or the static 
model, one other model in this family, corresponding t o the 
case F = introduced in I Miller fc Rouetl ll2002fl. ha s 



been studied in I Miller fc Rouetl 12006); Mill er et al 



iMiller fc Rouetl l|2010al ) as well as in liovce fc Sicar- 



We note, using Eq. @, that this class of models is ob- 
tained by taking the functional dependence of the 3D EdS 
expansion law, but allowing a freedom in the normalisation 
of the expansion rate to the matter density, i.e., we can ob- 
tain these models taking H 2 — K 2 87rGpo/3a 3 where k is a 
positive constant. Eq. Q then still holds, and thus a — e 2Tr , 
but instead of © we have, using Eq. @, that 



«Vffno/3. 



(6) 



In terms of 3D cosmological models the case K > 1, i.e., 
r > TEdS, can be considered to be a model in which there 
is not just "ordinary" matter but an additional pressure- 
less component of the energy density which does not clus- 
ter, e.g., a homogen eous scalar field w i th ap propriate equa- 
tion of state (see e.g. iFerreira fc Joyce] l|l998l ) and references 
therein). Equivalently one can consider that the k / 1 mod- 
els are obtained using the "correct" normalisation for the 
Hubble constant, H 2 — 8nGpo/3, but choosing the ID cou- 
pling freely instead of imposing the relation Eq. (|2]). It is 
important to note that, given that 3D Hubble law is im- 
posed "by hand" on the ID model, there is in principle no 
correct value for F. Our point of view here is simply to use 
r (or k) as a control parameter in trying to understand the 
role of expansion in non-linear structure formation. In other 
words we are seeking to improve our qualitative understand- 
ing only through the study of these toy models. 



3 STABLE CLUSTERING PREDICTIONS 

In the context of the problem of structure formation the 
family of models we study has the interest of being, like 
EdS models and static models in 3D, scale-free: the expan- 
sion introduces no additional time or length scales. This 
has as a consequence that, for cold particles with initial 
density fluctuations with a power law spectrum, one ex- 
pects to obtain, at sufficiently long times, temporal evo- 
lution which is "self-similar", i.e., the evolution in time of 
the spatial correlation functions is equivalent to a simple 
time dependent rescaling of the spatial variables. The rea- 
son why this occurs is that there is, in this case, only a 
single physically relevant length scale in the problem: the 
scale at which fluctuations are of order one, of which the 
evolution is predicted by linear theorjQ. Further if the clus- 
tering in the strongly non-linear regime does not depend 
on this scale one would expect to obtain a scale-free be- 
haviour of spatial correlations in this regime. Such power 

1 In the discrete problem there at least one additional length scale 
— the initial interparticle spacing — but physical (continuum 



law behaviour has been ob s erved in both expa nding (see, 
e.g., lEfstathiou et all (|l988l ); ISmith et all (120031) and refer- 
ences th erein) and static 3D simulations ( Baertschiger et ail 
l2007al lbl.l 20081 ) for a range of power law initial conditions. 3D 
simulations, however, can show such behaviour over a very 
limited spatial range, making it difficult to establish if it is 
associated with a true scale invariance of the non- linear clus- 
teri ng. The study of ID models in the family we are consider- 
ing (IMiller fc Rouetl2006l:lMiller et al.ll2007l ; IMiller fc Rouetl 
l2010al ; |Jovce fc Sicardll201ll l. with much greater spatial res- 
olution and numerical accuracy, show very convincing ev- 
idence that these models do indeed give rise to power-law 
clustering indicative of truly scale-invariant clustering. 

A central problem of non-linear structure formation is 
then that of understanding how the exponent (or possibly 
exponents) characterizing the non-linear cl ustering i s dete r- 
mined. In 3D EdS cosmology Peebles (see IPeebles! l|l98Ch ) 
proposed many years ago a simple explanation of how such 
scale invariant clustering might arise: if highly non-linear 
structures (of all different sizes) are supposed to evolve es- 
sentially independently after their formation, they will viri- 
alize and remain stable in physical coordinates. From this 
one may derive an exponent which depends only on the 
exponent in the power spectrum of initial density fluctua- 
tions. We derive he re now, more rigorousl y than the heuristic 
derivation given in Ijovce fc Sicardl l|201ll ) , the prediction of 
such a stable clustering hypothesis for the ID models we 
study. To do so we simply derive and then analyse the equa- 
tions of motion of a finite sub-system. 



3.1 Equations of motion of a finite subsystem 

In the simulations which we consider here, we start, as in 3D 
cosmological simulations, from a particle configuration gen- 
erated by applying small displacements to an infinite perfect 
lattice configuration. At any moment in time such a config- 
uration can be fully specified by giving the displacements Ui 
of particles from the sites of a perfect lattice. The position 
of the i-th particle with respect to some arbitrary origin, is 
Xi — iX + Ui, where A = rig . In a particle labelling with 
Xi ^ Xi-i for all i, the displacements are uniquely defined, 
and such that no two particles cross when t hey are applied 
to th e lattice. In this case it has been shown (|Gabrielli et all 
2009) that the force on the particle i, given by the sum on 
the right side of Eq. ((3]), with an appropriate specification 
of the regularisation, has the simple exact expression 



g sga(xi 



- Xj) = 2gn (ui — {u}) 



limit) results should not depend on it. 



where (u) is the average particle displacement. This result 
has been shown for an infinite system wi th statistically 
transl ationally invariant displacements in iGabrielli et alj 
( 2009), with the sole assumption of decay of correlations of 
the displacements at asymptoti cally large scales . The s ame 
result has been shown to hold in I Miller fc Rouetl (|2010bh for 
the case of an infinite periodic system, using an appropriate 
regularisation of an Ewald summation. In this case, as we 
will now show explicitly, the expression (O can be written 
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a.- 

j 

-9^2sgn(xi - xj) = -g[N L (i) - Nn(i)] + 2gn [xi - (x) p ] 

(8) 

where Nl(i) (Nr(i)) is the number of particles in the pe- 
riodic cell to the left (right) of particle i, and (x) p is the 
average particle position, i.e., centre of mass of the particles 
in the cel0. 

Let us consider now any finite subsystem S of the infi- 
nite system defined simply by the N particles, with i — 1...N 
say, in some interval [xo, xo + L] at a given instant in time, 
Substituting u% = Xi — iX in 0, and noting that the number 
of particles in the interval to the left (right) of the ith par- 
ticle is given by Nl(i) = i — 1 (Nr(i) = N — i), we obtain 
that the force on the particle i may be written 

F(xi) = -g[N L (i)-N R (i)]+2gno[xi-{x)s]-2gno[(u)-(u)s] 

(9) 

where (x)s is simply the position of the centre of mass of 
the subsystem, and (u)s the average displacement of the 
particles in the subsystem. If the subsystem chosen is the 
periodic cell of an infinite periodic system, we have (u) = 
{u)s and thus obtain ©. For any generic subsystem S, on 
the other hand, the last term on the right in © is the same 
for all particles in the subsystem. Defining yi = Xi~ {x)s we 
thus obtain the equation of motion for any particle in the 
subsystem, relative to the CM of the subsystem, as 

^ + r ^ = -9 E s g a(y i - Vj ) + 2gn y i (10) 

where the sum extends over the particles in the subsystem. 
This remains valid for as long as no other particle outside 
the system crosses the particles at its extremities. Eq. (|10p 
are thus exactly the equations of motion of a subsystem for 
as long as it is isolated in this sense. 

3.2 Definition of "quasi-physical" coordinates 

For r we can transform to the coordinates 



in which Eq. (|1U[I becomes 

^p=-9 2^ sgn(r 1 -r J ) + ^-[l + T ^]r l . (12) 

These coordinates have been chosen in order to remove the 
damping term in the equations of motion, and all depen- 
dence on the expansion now appears only in a time depen- 
dent force term, of which the amplitude decays as 1/P. At 
sufficiently long times, therefore, an isolated subsystem of 
the infinite ID "expanding" model has, to an arbitarily good 
approximation, the equations of motion of a strictly isolated 
non-expanding ID self-gravitating system. We will refer to 
the coordinates in |TT} as quasi-physical coordinates: they 
are the closest analogy we have to what are called physical 
coordinates in the three dimensional problem (i.e. f — a(t)x, 

2 Contrary to what is stated in lMiller fc R ouct (20l53), the ex- 
pressions for the force in the two papers are in agreement. 



t). In these coordinates the 3D cosmological equations of 
motion are equivalent to those for an infinite Newtonian 
self-gravitating system expanding about a centre in an in- 
finite (static) space. Correspondingly, any subsystem (e.g. 
galaxy) taken far from other mass in an expanding universe 
will evolve in th ese coordinates like an isolat ed Newtonian 
system (see e.g. IJovce fc Svlos Labinil (|2012T l for a deriva- 
tion). The additional term in 1/P in the ID model, which 
we cannot get rid of by a coordinate transformation, appears 
because we have imposed in the derivation of the model 
the 3D expansion law rather than using the analogy of this 
law in ID (which gives a collapse in a finite time, qualita- 
tively different to that of the 3D expansion which we wish to 
model). It is important to note that, re-expressed in terms of 
the three dimensional scale factor a(t), the transformation 
in (|11[) from the "comoving" coordinates yi to the "quasi- 
physical" coordinates in is r% = e 2TT ^yi = a 1 ^ 3 yi, and 
not Ti = ayi as one might naively expecfl 

3.3 Scaling of a finite virialized subsystem 

The second term on the right hand side of Eq. (|10[) may 
be rewritten simply as — h l]gnoyi The first term is of 
order gN ~ gnsL where ns is average comoving density 
of the subsystem. Thus if we consider any isolated subsys- 
tem which is already highly non-linear (with ns 3> no) we 
are necessarily in the regime where the second term may 
be neglected. In this case the isolated subsystem will thus 
behave, in quasi-physical coordinates, like a truly isolated 
ID self-gravitating system. Just as in 3D, such a system, for 
large N, in the collisionless limit, attains on a few dynamical 
time scales (~ 1/ \/gns) a s t ationa ry virialized state (see e.g. 
Ijovce fc Worrakitpoonponl (|201 lh and references therein), 
i.e., a system which becomes time independent in the quasi- 
physical coordinates r*. More generally, even if not strictly 
virialized, it rapidly attains a well defined average size about 
which virial oscillations persist. 

3.4 Linear amplification and self-similarity 

In order to derive the exponent in what follows we will also 
use the self-similar behaviour of the evolution, which, for the 
two point correlation function £(x, r) means that 

in the spatial range where the self-similarity is observed Q- 
The function R s (t) is derived from linear theory, which de- 
scribes a self-similar evolution from power law initial con- 
ditions once the growing mode dominates. In the ID model 
the linear theory growth of perturbations can be easily de- 
rived by using the expression for the force Eq. (0 in terms 

3 We note that lYano &: Go uda (1998) give an incorrect general- 
ization of the stable clustering hypothesis to the ID EdS model 
assuming that stability will be attained in the coordinates ax. As 
a result they arrive at the incorrect conclusion from their numer- 
ical study that stable clustering is not observed in the model. 

4 We follow the standard notation here in which £(x) de- 
notes the reduced two point correlation function, i.e., = 
n n 2 (n( x\)n(x\ - \- x)) — 1 where n(x) is the loca l particle density. 
See e.g. iPeeblesI dl98(t) or lGabrielli et al.l i2005h . 
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Figure 1. Two point correlation function measured in an evolved 
simulation, with k = 1 and n = 2. The values of x m i n and x m ax 
define the range in which there is (by eye) an apparent power-law 
behaviour. 



of the displacements 114. Setting, without loss of generality, 
the mean displacement (it) to zero, we have simply 

d 2 Ui diii 

__ +r _ =2^. (14) 
which is characterized by the growing mode 



u l (xe 2aTT =a a where a = ~[-l + Jl + ^] . (15) 

For k = 1 we thus recover the same growth law as in 
the standard 3D EdS model. Assuming a power law power 
spectrum of density fluctuations P(k) oc k n (and thus 
£ oc l/x n+1 at large x) the evolution is self-similar with 



Rs 



e i+„ 



0,1+n 



(16) 



3.5 



Exponent of non-linear two point correlation 
function 



Shown in Fig. [T] is a typical two point correlation £ (x) ob- 
tained from the evolution of power law initial conditions in 
the class of models we study. Between a lower cut-off scale 
Xmin, and an upper cut-off scale x„ lax a simple power law 
behaviour of the correlation function is observed £ ~ x~ J . 
The exponent 7 may be written as 

ln(£(a: m ax))-ln(£(3; m i n )) 

In (r ) - In (x ■ 1 ' ' 

The temporal evolution of the correlation function is 
found to be self-similar, down to the scale x m i n , below which 
it is broken. The scale a; m ax is in the range in which self- 
similarity applies, and thus scales in proportion to R s , with 
£ (£max) a fixed amplitude corresponding to the transition 
between the linear and non-linear regime. 



If we make now the hypothesis that the structures giv- 
ing rise to the correlation measured at the lower cut-off are 
stable in quasi-physical coordinates, i.e., these structures be- 
have like isolated approximately virialized structures, we 
have that x m i„ ~ a~ 1,/3 . Further, since for £ 3> 1, £(x) 
mesures simply the average density at distance x from an 
occupied point normalized to the mean density, we have that 



£ (^min) 



,V3 



It follows that the exponent tends asymp- 



totically (at sufficiently long times) to 



7sc (n, k) 



2k (n + 1) 



(18) 



«(2n- 1) + 3V^ 2 +24 

This result is identical to that derived in Ijovce fc Sicardl 
l|201ll ) using more heuristic arguments based on the scaling 
of the energjQ It tells us how the exponent, in the sta- 
ble clustering hypothesis, depends on the initial conditions, 
characterized by the exponent n of the power spectrum of 
density fluctuations, and on the expansion rate parametrized 
by the constant k, which we recall is 



K = Ts/3/gno = f ^ " 
\87rGpo 



2 3\ 1/2 



(19) 



Note that as k — > we obtain 7 SC = 0, i.e. a correlation 
function which is predicted to be flat over an arbitrarily 
large distance. Indeed, in this static limit, the stable clus- 
tering hypothesis corresponds to stability also in comoving 
coordinates, and therefore to a constant x m in and £ (xmin). 
Clearly, however, we do not expect in this case that the for- 
mation of non-linear structure at a given scale will leave ap- 
proximately unperturbed the structures formed previously 
at smaller scales. 



NUMERICAL SIMULATIONS: METHODS 
AND RESULTS 



4.1 Method and calibration 



The numerical study in IJovce fc Sicardl l|201ll ). as most pre- 
vious studies of the three models belonging to the class we 
are considering, exploits fully the very attractive property 
of these models that they admit "exact" numerical integra- 
tion: the force on particles, whether given in the form ((7} or 
JS), is such that the motion (with or without damping) ad- 
mits an analytical solution other than at particle crossings. 
An event driven algorithm can then be employed, in which 
the determination of the time of the next "collision", and 
which pair of particles it involves, requ ires only the solution 
of algebraic equations. As described bv lNoullez et al.l (|2003l ) 
the integration can be sped up optimally using a "heap" 
structure. 

For the previously considered cases, corresponding to 
k — 0, l,\/3, the algebraic equations involved are, respec- 
tively, quadratic, quintic and cubic. For a generic k, however, 
the equations are not polynomial, and the numerical cost of 
solving them increases and makes the code, which is already 
expensive numerically, even more time consuming. As we 
wish to study fully a broad range of k and initial conditions, 



5 In I Joyce fc Sicardl fcullTl . as in various other papers on these 
models, units were used in which 2gno = 1, and thus Y = ft/\/6. 

6 The equation for the crossing times is of the form 1 + Az& + 
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we have chosen here instead to use a particle-mesh (PM) 
code of which the numerical speed can be greatly enhanced 
by the use of FFT techniques and parallelizatioqj. In order 
to be sure that we are not, as a result, losing the advantages 
of the precision of the integration accessible in these ID 
models, we calibrate, as discussed below, our choice of the 
new discreteness parameters introduced by the PM code by 
comparing our results with those obtained, for k — 0, 1, v3, 
using the "exact" code. 

Our PM code defines a regular grid on the simulation 
box, on which a density field n(x) is defined at each time step 
by a "clouds in cell" interpolation of the particle positions. 
The modified ID Poisson equation 



d 2 u 



= 2g[n{x) - n ] 



(20) 



is then solved for the potential in Fourier space using an 
FFT, using the solver fftw3 0. An inverse FFT then deter- 
mines the force on the grid, which is then interpolated back 
onto the particle positions. The latter are then advanced us- 
ing a leapfrog algorithm, with a simple adaptative timestep: 
it is chosen at each time by imposing that the mean distance 
travelled by particles during a time step be equal to a chosen 
fraction r\ of the PM grid. 

Compared to the exact code, our PM code therefore in- 
troduces resolution effects controlled by two dimensionless 
parameters: the size of the PM grid compared to the lattice 
spacing in the initial conditions, which we denote by e, and 
the parameter r\ controlling the size of the time steps. We 
have made an extensive study of the effects introduced by 
this finite resolution of the PM code, and in particular on the 
determination of the two point correlation function which is 
the quantity which interests us here. Shown in Fig. [5] is, for 
example, a comparison of the results for runs from identical 
initial conditions, for the case n = 2 and K = 1, obtained 
from the exact code and the PM code, for two different val- 
ues of e, and r\ sufficiently small (of order unity) so that 
convergence is obtained with respect to it. We find, very rea- 
sonably, that excellent agreement is obtained provided e is 
chosen smaller than x m i n , while a larger e than this minimal 
scale leads to a visible deviation of the correlation function 
from its correct value at scales below e. We have chosen our 
numerical parameters in the simulations reported below in 
order to obtain such an indiscernible difference in the cor- 
relation function between our code and the exact code, for 
the cases k = 0,1, V3. All the results presented here are 
for systems with N = 10 5 particles (in periodic boundary 
conditions). 

For initial conditions with power law power spectra 
P(k) ~ k n (at small k, a cut-off at large k is always im- 
plicitly assumed), we expect hierarchical structur e forma- 
tion in d dimensions for n in the range —d < n ^ 4 ij Peebles! 
1980): for n ^ — d the variance of mass fluctuations diverges 
at large scales, while for n > 4 linear theory breaks down 
(and the dynamics at large scales will be determined by 
the power at the ultra-violet scale). We thus consider (as in 



Bz (p+l)/2 = o w here z = e rr , /3 = Jl + and A and B are 



constants 
7 



See, e.g., |http://astro.uchicago.edu/^ andrey/Talks/PM/p m.pdf| 
and references therein. 
8 See http://www.fftw.org/ 
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Figure 2. Two point correlation functions obtained by evolving 
from the same initial condition (with n = 2) and re = 1, for three 
different cases: using the "exact" code (solid line) and using our 
PM code with the two indicated values of the PM grid spacing. 



iJovce fc Sicardl (^OllD ) the cases n = 0, 2, 4. Point processes 
with such power spectra are easy to set up: n — is obtained 
by randomly distributing points in the interval (giving a 
power spectrum P(k) — 1/no), the case n = 2 by applying 
small random and uncorrelated displacements to a regular 
lattice, and n = 4 in the same way but with the additional 
constraint that nei ghbouring pairs have equ al and opposite 
displacements (see iGabrielli fc Joyed 1)20081 ) and references 
therein) . 



4.2 Results 

Our results are summarized in Fig. [3] Each point on the 
plot corresponds to a single simulation with given n and 
k, where the latter varies in the range — 2.5. The figure 
shows then the value of the exponent 7 of the power law re- 
gion of the correlation function. In each case the given value 
is determined using a linear regression in the range between 
i m ,„ and Xmax, the latter being determined themselves by 
eye (cf. Fig. [TJ. The uncertainty in the fitted exponent as- 
sociated with the (by eye) estimation of the end-points for 
the fit is very small (of order a few percent at most) in 
most parts of the parameter space. However, as we will now 
discuss and explain, at lower values of k and n, where the 
extent of the power law region is much smaller, the resul- 
tant uncertainty becomes more significant. The continuous 
curves plotted correspond to the predictions of the stable 
clustering hypothesis, Eq. {IS) . 

Our results show excellent agreement with the stable 
clustering prediction for a very large part of the explored 
parameter space of initial conditions and cosmology. This 
region can be well characterized simply by the condition 
7 sc (n, k) > 0.2, i.e., good agreement with stable cluster- 
ing is observed simply in the region where the predicted 
exponent be larger than some critical value. In the rest of 
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Figure 3. Exponents of non-linear power law clustering 
7 measured in simulations with different initial conditions 
(parametrized by n) and expansion rates (parametrized by k). 
The solid lines are the predictions of stable clustering, Eq. JTS). 



Figure 4. Schematic representation of our results: the space of 
initial conditions (n) and cosmology (k) breaks into two parts, 
with a boundary defined approximately by 7 sc (n,re) = 0.15. 



the parameter space (i.e. for y sc (n,K,) < 0.2) there is clear 
disagreement with the stable clustering prediction, and the 
measured exponents lie in a very narrow range, between 0.2 
and 0.15, or a little smaller. It is precisely in this part of the 
parameter space, however, that there is also a considerable 
scatter about the curves of the measured exponents. This 
scatter is in fact just a measure of the uncertainty in the de- 
termination of these exponents, which is most difficult (for 
reasons we explain in detail below) in the region where the 
exponents become small. Thus, our results are quite con- 
sistent with the hypothesis that the exponent in the region 
7 sc (n, k) < 0.2 is "universal", in the sense that it is inde- 
pendent of initial conditions and cosmology. A schematic 
representation of this result is shown in Fig. 0] the parame- 
ter space of initial conditions and cosmology breaks up into 
a "stable clustering" region, and a region of "universality" . 
The contour defining the two disjoint regions is simply given 
1sc(n, k) = 7o- Rather than a line, there may of course be a 
small region in parameter space in which the exponent varies 
in a continuous manner from its stable clustering value to a 
limiting universal value of 70 ~ 0.15 0. 



5 DISCUSSION AND CONCLUSIONS 
5.1 Interpretation of ID results 

We have found numerically that the stable clustering hy- 
pothesis predicts the exponents of non-linear clustering very 
accurately at sufficiently large n and k, i.e., if the spectrum 
of initially fluctuations is "sufficiently blue" and the expan- 
sion rate is sufficiently fast. More precisely it appears that 
the criterion for it to work is that 7 SC , which is a mono- 
tonically increasing function of both n and k, be larger than 
some critical value. It is not difficult to give a simple physical 
explanation for this behaviour, as we now explain. 

Let us consider two overdense regions, initially (in the 
linear regime) of comoving size L\ and L% > L?. Self- 
similarity of their evolution implies that there is simply a 
time delay Ari2 between every stage of their evolution given 
by R s (Ar 12 ) = (L%/Ll). Using Eq. t\W$ we obtain 



Ari2 



1 + 71 

4Fa 



(21) 



Let us suppose that the fluctuation 1 goes non-linear and 
approximately virializes at some given time t\ , when its size 
is L\. When fluctuation 2 virializes, at time t% = T\ + Ari2, 
its size is L% = (L\l L\~)L\ . If in the interval between r" and 
t%, structure 1 behaves to a very good approximation as 
an isolated structure, its size in comoving coordinates will 
decrease by a factor e _2rAT//3 . It follows that the relative 
size of the two structures at the time when 2 virializes is 



9 In the language of statistical physics our result could be de- 
scribed as an evidence for an "out of equilibrium phase transi- 
tion" : as control parameters for the system are changed the (out 
of equilibrium) behaviour of the system changes in a discontinu- 
ous manner. With the precision of our current results we cannot 
determine with confidence whether there is really such a transi- 
tion or rather a continuous "crossover" between the two "phases" . 



(22) 



where a is the growing mode exponent given in ((15} . Quite 
simply, the faster the expansion rate, or the larger is n, the 
more a given scale which has entered the non-linear regime 
can "shrink" relative to a larger scale in the time before 
the latter goes non- linear. Further we note, using ([15}, that 



© 0000 RAS, MNRAS 000, 000-000 



8 D. Benhaiem, M. Joyce and F. Sicard 



Eq. 



can be written 



(23) 



i.e., the dependence on (n, k) is in fact completely deter- 
mined through the exponent 7 SC itself. 

The predicted stable clustering exponent is therefore di- 
rectly related to a physical quantity which we would expect 
naturally to control the validity of the assumption made 
in deriving it: the more "concentrated" are the pre-existing 
virialized substructures inside a larger structure when it col- 
lapses, the better should become the approximation that this 
structure will behave simply as a collection of sub-structures 
which are not disrupted by the virialization and subsequent 
evolution of the larger structure. Indeed, in the limit that 7 SC 
approaches unity, any structure which collapses and virial- 
izes will see the substructures it contains essentially as point 
particles. As 7 SC decreases, on the other hand, we expect 
that the interaction between structures can lead to their 
disruption, and in particular that "mergers" of substruc- 
tures become much more probable. Given that the value of 
the exponent 7 SC characterizes precisely the relative "con- 
densation" of scales, it is very natural that a critical value 
of the exponent should characterize the breakdown of stable 
clustering. 

What our results do not allow us to conclude, as we have 
discussed, is what precisely happens in the regime where 
stable clustering breaks down: our numerical results do not 
allow us to clearly distinguish between the possibility of an 
abrupt (discontinuous) transition to a region in which there 
is a truly universal exponent, or a smoother transition to a 
region in which the exponent depends only very weakly on 
the relevant parameters (n, n). Further we do not (currently) 
have a model to explain the exponent (or narrow range of 
exponents) observed in this part of the parameter space. 

Extending the analysis given just above to determine 
the scalings of relative sizes of structures, it is simple to un- 
derstand why it is precisely the region of small 7 SC in which 
the numerical results for the exponent of the power law in 
the correlation function are most noisy. The degree of preci- 
sion in the measured exponent is essentially just a function 
of the range of scale over which the power-law behaviour ex- 
tends, i.e., it depends on the ratio x m i n /x max accessible in a 
numerical simulation. The scale Xmax corresponds approxi- 
mately at any time to the size of the largest (approximately) 
virialized scale. If we call x% the comoving size of the first 
structure which virializes in the simulation, at the time r, at 
which it virializes. Self-similarity and stable clustering then 
imply the temporal evolutions 



x min ^ ^_2r(T-T )/3 



Rs(T) 

Rs(r ) 



(24) 

We can infer, using Eq. (|2ip. that, the range over which 
power law scaling is expected for r > to at the end of the 
simulation can be expressed as 



In 



1 — 7s, 



•hi 



Li 



1 — 7s, 



In 



(25) 



where L° and L® are, respectively, the initial comoving size 
of the first and last scale to virialize during the evolution, 



and Nf (JVj) are the number of particles they contain, re- 
spectively. Thus, while for 7,. mac closer to unity this range 
of scales is large, as 7 SC decreases it contracts and is simply 

L° 

of order the ratio jfr for 7 SC w 70. 

For our simulations, with TV = 10 5 particles, taking 
Ni ~ 10, or a little larger, and Nf ~ 10 3 — 10 4 , we can, for 
small 7sc, access power law clustering over approximately 
two orders of magnitude, which leads in practice to the ob- 
served order of the uncertainly in the exponent (~ 0.05). 



5.2 Comparison with 3D simulations: correlation 
analysis 

Let us compare our results to those in three dimensions. 
Numerical simulations have been done by various groups to 
test the stable clustering hypothesis, via the study of the 
EdS model starting f r om po w er-law initial conditions (see , 



Efstathiou et al 



(1988); Padmanabhan et al 



Colombi et alJ (11999): I Jain feBertschingerl (|l996 



1996) 



1998) 



Ma fc Frvl (|2000bl) ; ISmith et all (120031 )). In this case, corre- 
sponding to k — 1 in our notation, the predicted exponent, 
in three dimensions, of the two point correlation function is 
73 = 3(3 + n)/(5 + n). The range of power spectra which 
has been probed numerically is —3 < n ^ 1. Further, more 
generally, the stable clustering hypothesis leads to scaling 
relations for higher order co rrelation functions, which have 
been t ested by some authors (|Colombi et alJl996tlMa fc Frvl 
l2000bl ). 

Until the m ost recent and extensive study of 
ISmith et al.1 (|2003l ). previous works had concluded that the 
stable clustering prediction worked quite well, with, in some 
cases, indications of small deviations close to the limits of 
numerical resolution. This last study, on the other hand, 
reports measurements of the exponent 7 which are clearly 
discrepant with the stable clustering prediction: for n = —2, 
7 = 0.77 rather than 73 = 1; for n — —1.5, 7 = 0.91 rather 
than 73 = 1.29; for n — — 1, 7 = 1.26 rather than 73 = 1.5; 
for n — 0, 7 = 1.49 rather than 73 = 1.8. 

These results are clearly qualitatively different from 
those we have found in our ID models, for any n. While 
the measured 7 increases as n increases, just as in one di- 
mension, the discrepancy with stable clustering in three di- 
mensions is associated with an exponent 7 which is smaller 
than the predicted stable clustering exponent. Further there 
is no indication — in this limited range of n probed — that 
the discrepancy with respect to stable clustering is increas- 
ing, nor that the associated discrepancy manifests itself as 
an (exact or approximate) universality of the measured ex- 
ponent. 

These differences between ID and 3D results could, of 
course, be simply a reflection of the difference of the grav- 
itational clustering in the two ca ses. One difference wh ich 
may be essential, as discussed in Ijovce fc Sicard] l|201ll ). is 
that in one dimension there are no tidal forces, and a virial- 
ized structure may only be perturbed by another structure 
actually crossing it (which is, on the other hand, much more 
likely than in three dimensions). However, one of the main 
motivations for this study is precisely that the reliability of 
results from 3D simulations is not clear: the combination of 
the fact that the range of comoving scales probed grows in 
proportion to N 1 ^ 3 , and the necessity to introduce a (rel- 
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atively large) smoothing parameter to regularize the small 
scale divergence in the 3D gravitational force, mean that the 
exponents (in existing studies) are measured over at most 
a little more than one order of magnitude. Further there 
are consideable discrepanc ies between some of the studies. 
iJain fc Bertschingerl (|l998l ) find, for example, good agree- 
ment with stable clustering for th e case n = — 2 with simu- 
lations of a similar size to those of ISmith et al. | l|2003lf °l 

Let us focus a little more on these differences between 
ID and 3D simulations linked to the small scale smoothing 
of the force. In ID, as discussed, the equations of motion can 
be integrated exactly in absence of smoothing, so that there 
is in practice no lower limit on spatial resolution other than 
that arising from the finite particle density, of order the scale 
Xmin in the correlation function. We note that, following the 
analysis above, we infer that, if stable clustering holds, we 
have 



In 



1 — 7sc 



hi 



(26) 



in one dimension, while the same arguments applied to three 
dimensions give 



In 



7s, 



hi 



7sc 



(27) 



The scale Xi, the comoving size of the first structure which 
virializes, is of order the initial interparticle spacing A. In 
one dimension, even when we use, as here, a PM code, the 
smoothing can, without excessive numerical cost, be chosen 
so that it is at all times considerably smaller than the scale 
Xmin, so that we can be absolutely confident that it plays no 
role in the evolution. In three dimensions, on the contrary, 
cosmological simulations of a reasonable size (i.e., so that 

L°, 

the range of scales jjj which go non-linear cover at least a 

decade or two) are typically performed with a smoothing of 
between one tenth and one hundredth of the interparticle 
spacing A. Thus 3D simulations cannot in practice simu- 
late systems manifesting stable clustering only over a very 
limited range of length scales, while remaining in the regime 
where e <C x m i n . While the use of a cut-off e ^ x m in does not 
necessarily imply that the clustering above the scale e is not 
accurately reproduced, it is quite possible that this could be 
the case if great care is not taken in t he choice of numerical 
parameters (see iKnebe et al.l (|2000h : iJovce fc Svlos Labinil 
20121 )). In this respect even when a relatively large smooth- 
ing is employed, two body relaxation may play a role in dis- 
rupting the desired collisionless evolution of iV-body simula- 
tions, and such effects c an potentially lead to a "pollution" 
of larger scales than e (|Knebe et al. l200d) . In the ID sys- 
tem, in contrast, there is in fact no two body collisionality 



10 This notable difference is attributed by t he latter, on the ba- 
sis of a qualitative analysis, to the fact that I Jain fc Bertschingerl 
( 1998) evolve to a time at which the amplitude of the power at the 
scale of the box is higher, leading possibly to poor modelling of the 
coupling to the "missing" long wavelength modes. Another no- 
table difference in the simulations is, howeve r, that a considerably 
smalle r smoothing length is employed by Ijain fc BertschingeJ 
(1998). As discussed below, a larger smoothing could be respon- 
sible for suppressing power on smaller scales, leading potentially 
to an effective exponent which is lowered. 



analagous to that in three dimensions, and collisional relax- 
ation (of which the mechanism is not yet fully understood) 
is extremely sl ow compared to that induced by 3 D two body 
relaxation (see IJovce fc WorrakitpoonporJ (|2010h and refer- 
ences therein). Simple estimates show that such relaxation 
should play absolutely no role in our ID simulations, while 
the same is not true in typical 3D simulations. Our study 
thus suggests that particular attention should be paid to 
this point in the analysis of 3D simulations. 

We note also that simulations of clustering in a 3D 
universe have been perfor med for the case of a static uni- 
verse, for the cases n = (iBottaccio et al.ll2002h and n = 2 



|Baertschiger et all l2007al lbl l200&t ). The results are in this 



case q ualitatively very consis tent with those observed here 
and in IJovce fc Sica rd (2011) for the static (k = 0) limit. 
Just as in the ID simulations self-similarity is observed, and 
further a correlation function which appears to be the inde- 
pendent of the initial conditions, with a very shallow expo- 
nent in the inner part, 7 « 0.3. 

5.3 Comparison with 3D simulations: halos and 
universality 

Finally let us discuss the question of "universality" . A strik- 
ing feature of our ID results is that the measured exponent 
characterising clustering appears to become independent of 
initial conditions and cosmology when the stable clustering 
approximation breaks down — more specifically when the 
exponent n in the initial power spectrum is sufficiently small 
for any given expansion rate paramet er k. In three d imen- 
sions , following notably the work of iNavarro et al. I (Il996l . 
1997), an apparent "universality" of gravitational clustering 
which manifests itself as an independence of halo profiles 
of initial conditions and cosmology. The range over which 
such a universality might apply is not, how ever, clear: stud- 
ies such as that of iKnollmann et al.l l|2008l ) show that, even 
in the EdS cosmology, the inner exponent of halos may vary 
sensibly for power law initial conditions as a function of n. 

Does the apparent universality we observe in ID cor- 
respond to a u n iversa lity in halo profiles? As discussed in 
IJovce fc Sicardl i|201ll ). one can identify and extract halos 
by prescriptions analogous to those used in three dimen- 
sions. However, the objects thus extracted are simply so 
"clumpy", at any scale significantly above the scale x m i n , 
that it is not possible to approximate the run of density by 
a smooth profile: indeed it can be verified directly that the 
power law exponent (over several decades) in the correla- 
tion function corresponds to a truly scale invariant fractal 
distribution of the matter in the corresponding range. If, 
however, the extracted halos are smoothed (by hand) over 
a scale much larger than x m in, the exponent characterising 
the decay of density about the centre of the halo, should be 
approximately the same as that measured in the correlation 
function. This is true because the latter measures, by defi- 
nition, the mean density as a function of distance about any 
point, which in a simple fractal type distribution is (mod- 
ulo fluctuations) the same about any point. Thus, in the ID 
model, the universality we observe here is indeed expected 
to be associated to a universality of inner exponents of halos. 

Extrapolated to 3D, our results thus suggest that the 
relative smoothness of halo profiles is an artefact of insuffi- 
cient spatial resolution, and that the minimal scale of sub- 
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structure is limited only by such resolution. In this case, for 
power law initial conditions, the observed power law slope 
of the correlation function should in fact coincide with the 
inner exponent of the measured halo profiles. We note that 
this is completely different to what is predicted usually in 
thre e dimensio n s on th e ba sis of a smooth h alo m odel (see 
e.g. iMa fc Frvl |2000bl lcj|ah; lYano fc Goudal |2000r i). Up to 
now 3D numerical studies have not, to our knowledge, ad- 
dressed this point directly (by comparing measured slopes 
in halos and two point correlation functi ons, in scale free 
cosmo logies). We note that the results of iKnollmann et al.l 
(2008) on inner halo profiles appear, however, to be very 
consistent with a behaviour completely analogous to what 
we have seen in our ID models: in their study of power law 
initial conditions in an EdS cosmology, for n varying in a 
range between —0.5 and —2.75, inner slopes for halo pro- 
files are estimated using different fitting procedures. For the 
simplest such procedure it is found (see Fig. 2) that the mea- 
sured exponent varies from a value in good agreement with 
stable clustering for the largest value of n, but appears to 
deviate, as n decreases, towards an exponent which is larger 
than predicted by stable clustering, tending rapidly towards 
an asymptotic exp onent a little larger t han unity, and con- 
sistent with that of lNavarro et al.l l| 19971 ). This suggests that 
in the 3D EdS model, the universality of halos would trans- 
late to an insensitivity to initial conditions below an n of 
order —2. Given that the range of effective exponents of ini- 
tial conditions of typical cosmological (e.g. ACDM) initial 
conditions falls in or close to this range, this would appear 
to be very consistent with the approximate universality of 
halo profile exponents observed in these kind of simulations. 

5.4 Future studies 

Our ID study thus motivates further careful numerical study 
of non- linear clustering in the pure matter (k = 1) EdS 
m odel in three dimen sion s, extending recen t stud ies such 
as ISmith et~aH (|2003h and IKnollmann et all (|2008l ). In the 
same way as has been done here, the study of this model can 
be extended to the full family of EdS models of which one 
limit is a static universe. The goal of such a study would be 
to see if, in this space of models, and for power law initial 
conditions, a similar structure is found to that we have seen 
in one dimension. Our ID study suggests that particular at- 
tention needs to be applied to controlling for the robustness 
of measure of exponents in the two point correlation function 
to small scale force smoothing, particularly when an expo- 
ne nt below the predic ted stable clustering value is found (as 
in ISmith et al.l l|2003l )). Further, study of the relation be- 
tween exponents measured from the two point correlation 
function and those measured from halos, and their variation 
as a function of n and k, should clarify whether the qualita- 
tive nature of clustering is indeed like that usually assumed 
in theoretical modelling (smooth halo models) or instead like 
that in the ID models (hierarchical fractal clustering). 

If the framework suggested by the ID models — of a 
"universal" region in the space of initial conditions of clus- 
tering, delimited by a "critical" value of the predicted sta- 
ble clustering exponent — is correct, and applies also in 
three dimensions, the theoretical problem remains of un- 
derstanding both the corresponding value of the "universal 
exponent" , and why this universality applies in the specific 



region where it does. Further study of the ID model may 
then help to throw further light on this. As the static model 
is in the relevant part of the parameter space, it may suffice 
to study this particular limit, i.e., the region in which there 
is universality is where the expansion of the universe has lit- 
tle effect on the non-linear c l usteri ng. We note that a very 
recent study in ISchulz et alj (|2012l ) shows that there is an 
apparent universality in the properties of the "halos" formed 
from the ID collapse of a single structure with a range of 
cold initial conditions, and argues that this universality may 
be related to tha t observed numerica lly in three dimensions. 
The ID halos in ISchulz et all i|2012l) , however, are smooth, 
very different to the non-linear structures we observe here, 
which are very clumpy even in the static limit of the ID 
model. As discussed above, however, in this limit the range 
of non-linear clustering explored by our simulations is still 
quite modest, and it may be t hat larger simulatio ns might 
establish a link to the study of ISchulz et al.l l|2012T ). 

We are indebted to B . Marcos for numerous discus- 
sions and suggestions. We also thank S. Colombi, B. Miller, 
J. Morand, J.-L. Rouet, F. Sylos Labini, P. Viot and T. 
Worrakitpoonpon for useful conversations. 
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